Project-oriented workflow

As you know, we have created the RStudio project for AGE exercises, with its top-level directory being ~/AGE2021/Exercises (you can find the Exercises.Rproj file there). Whenever you open a project in RStudio, it will automatically set the working directory to its top-level.

Now that you have this information, you can think of “I will stop using setwd() and always use paths relative to top-level directory”. Unfortunately, some packages are evil and set your working directory silently, which leads to unexpected and hard-to-debug errors (one example is knitr described in intro_to_Rmd.Rmd). Luckily, for a consistent navigation in your projects, here package was created.

Let’s see how it works. The main function is here::here():

library(here)
here()

## [1] "/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises"

here("E02-intro_to_R", "images")

## [1] "/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises/E02-intro_to_R/images"

It joins paths relative (E02-intro_to_R/images) to the project’s top-level directory with its absolute path (/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises). But what happens if we change the working directory?

# Just for educational purposes!
setwd("E01-intro")
cat("Current working directory:\n")
getwd()

Current working directory:
[1] "/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises/E01-intro"

cat("here():\n")
here()

here():
[1] "/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises"

cat("Some path:\n")
here("E02-intro_to_R", "images")
setwd(here())

Some path:
[1] "/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises/E02-intro_to_R/images"

Hmm, nothing happened after setwd(), right? This way, no hidden setwd() can surprise you 😎

How is it possible here always knows where is the top-level of project located? Well, here is smart: it scans your project directory structure and do some heuristics. The most basic one is searching for file ending with .Rproj. Wait! Is that one of the reasons why we have created our RStudio project? Yes, exactly 🤓

If you want to save space and typing, you can also define variables with here paths and use them to construct paths relative to them:

output_dir <- here("some/output/directory")
here(output_dir, "analysis_1")

## [1] "/data/persistent/jirinovo/bio-class-deb10/AGE2021/Exercises/some/output/directory/analysis_1"

You can find more information about project-oriented workflow here and about safe here paths here and here.

Put runtime information to the end of your script/Rmd file

It is a good practice to do so for the sake of reproducibility. This code outputs all important information about your R session (see the Cleanup section for actual rendered output):

• Warnings.
• Traceback of errors.
• Loaded packages and their versions, external libraries (BLAS etc.).

warnings()
traceback()
sessionInfo()

Instead of sessionInfo() you can utilize the sessioninfo package, which offers a prettier output:

sessioninfo::session_info()

Do not automatically reload your R session

Otherwise you can expect the unexpected.

Always start R with a blank state.

Because you are using R through RStudio Server, don’t worry when you close your browser - R session will remain running.

Prevent package namespace conflicts

When you load a package with library(), all exported objects from that packages are attached to the global namespace, while replacing existing objects. Let’s see an example (restart your R session before). Function named select() is very common, for example it is present in these two packages:

AnnotationDbi::select

## standardGeneric for "select" defined from package "AnnotationDbi"
## 
## function (x, keys, columns, keytype, ...) 
## standardGeneric("select")
## <bytecode: 0x55c77de90da0>
## <environment: 0x55c77debd4e8>
## Methods may be defined for arguments: x
## Use  showMethods("select")  for currently available ones.

dplyr::select

## function (.data, ...) 
## {
##     UseMethod("select")
## }
## <bytecode: 0x55c76da0df80>
## <environment: namespace:dplyr>

If we use package prefix to locate the select() function (AnnotationDbi:: or dplyr::), no problems can arise.

But see what happens when we load these two packages:

library(AnnotationDbi)
library(dplyr)
select

function (.data, ...) 
{
    UseMethod("select")
}
<bytecode: 0x5558f010cba8>
<environment: namespace:dplyr>

You see that select() from dplyr won, because this package was loaded after AnnotationDbi. This could lead to unexpected results: you think you are using a function from a package, but instead you are using a function with the same name, but from an another package, which was loaded last.

To prevent such confusions, we can use the conflicted package. It is a sort of guardian and dictator in one person: it throws error when you attempt to use a function with ambiguous name, and also allows you to select a preferred package for such conflicts. Let’s see it in action (restart your R session again):

library(conflicted)
library(AnnotationDbi)
library(dplyr)
select()

Error: [conflicted] `select` found in 2 packages.
Either pick the one you want with `::` 
* dplyr::select
* AnnotationDbi::select
Or declare a preference with `conflict_prefer()`
* conflict_prefer("select", "dplyr")
* conflict_prefer("select", "AnnotationDbi")
Run `rlang::last_error()` to see where the error occurred.

You see that select is conflicting and conflicted offers you two options how to resolve it:

• Use the package prefix:

dplyr::select

• Or use conflict_prefer() to select a preferred package:

conflict_prefer("select", "dplyr")
select

[conflicted] Will prefer dplyr::select over any other package

function (.data, ...) 
{
    UseMethod("select")
}
<bytecode: 0x56350f4713e8>
<environment: namespace:dplyr>

conflicted should be loaded as the first package.

Unfortunately, I encountered some RStudio freezes during the usage of conflicted. In that case refer to the section “What to do when RStudio is stuck” in E01-intro.

Alternatively, the following packages offer a clearer mechanism of package loading. They allow to import only specific objects from a package, similar to Python’s from package import function as fn

• box
• import

Use renv to capture your project dependencies

This is really not needed for AGE, but anyway it is a good practice to somehow export all packages needed to run your project, and renv is a great tool for that. It basically does two things:

• Creates a separate R library for your project. This way, all packages are installed into this library instead of the system one. So when some project needs a particular package version, you won’t mess up your system library.
• It searches your project files for library() and package:: usages and exports packages with exact versions to a file called renv.lock. Based on this file, everyone is able to install all packages needed for a project.

renv is somehow similar to virtual environments in conda 🙂

Use targets: a Make-like pipeline toolkit for R

The targets package is a Make-like pipeline toolkit for statistics and data science in R. With targets, you can maintain a reproducible workflow without repeating yourself. targets skips costly runtime for tasks that are already up to date, runs the necessary computation with implicit parallel computing, and abstracts files as R objects. A fully up-to-date targets pipeline is tangible evidence that the output aligns with the code and data, which substantiates trust in the results.

This is really out of scope of AGE, but targets is the way reproducible R analyses should be done.

traceback()

When error occurs, its message is usually not very informative. traceback() function prints a call stack of last error, helping you to resolve a location of the error.

Let’s define the group of functions, which call each other, starting from f():

f <- function(x) g(x)
g <- function(y) h(y)
h <- function(z) {
  z**2
}

Now we can try calling f(). It will fail for "2", obviously:

f(2)

## [1] 4

f("2")

## Error in z^2: non-numeric argument to binary operator

We can inspect the call stack:

traceback()

3: h(y) at <text>#2
2: g(x) at <text>#1
1: f("2")

Also try running f("2") in your prompt. RStudio will directly offer you to show traceback or rerun code with debug (see below).

Interactive debugger

Sometimes you need more information from a location of error. Interactive debugger allows you to pause execution of a function and interactively explore its state. It is invoked by calling the browser() function:

f2 <- function(x) g2(x)
g2 <- function(y) h2(y)
h2 <- function(z) {
  browser()
  z**2
}

f2("2")

Interactive debugger can be also invoked by RStudio, as is shown above in the traceback() section.

Unfortunately, browser() is not working well from within code chunks, so rather run debugged function from the R command line (interactive session).

recover()

Another way to activate browser() is to use options(error = recover). Now when you get an error, you will get an interactive prompt that displays the traceback and gives you the ability to interactively debug inside any of the frames:

options(error = recover)
f("2")

To restore the default on-error behaviour, use dump.frames:

options(error = dump.frames)

Vectorized functions

Some of the examples are mean(), sum(), cumsum() or vectorized if ifelse():

my_vector <- 1:5
mean(my_vector)

## [1] 3

sum(my_vector)

## [1] 15

cumsum(my_vector)

## [1]  1  3  6 10 15

ifelse(my_vector < 2, "< 2", ">= 2")

## [1] "< 2"  ">= 2" ">= 2" ">= 2" ">= 2"

Apply-family functions

These functions take data (vector, list, dataframe, …) as an input and apply a function on each of its element. Let’s introduce the most basic ones:

sapply() apply a function to each element of a vector or list and returns a vector:

sapply(my_vector, function(x) {
  x ** 2
})

## [1]  1  4  9 16 25

lapply() is similar, but returns a list instead of vector. It can be also used for dataframes, where a function is applied to each of dataframe’s column (dataframe is basically a list of column-lists):

lapply(my_vector, function(x) {
  x ** 2
})

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 4
## 
## [[3]]
## [1] 9
## 
## [[4]]
## [1] 16
## 
## [[5]]
## [1] 25

lapply(mtcars, head)

## $mpg
## [1] 21.0 21.0 22.8 21.4 18.7 18.1
## 
## $cyl
## [1] 6 6 4 6 8 6
## 
## $disp
## [1] 160 160 108 258 360 225
## 
## $hp
## [1] 110 110  93 110 175 105
## 
## $drat
## [1] 3.90 3.90 3.85 3.08 3.15 2.76
## 
## $wt
## [1] 2.620 2.875 2.320 3.215 3.440 3.460
## 
## $qsec
## [1] 16.46 17.02 18.61 19.44 17.02 20.22
## 
## $vs
## [1] 0 0 1 1 0 1
## 
## $am
## [1] 1 1 1 0 0 0
## 
## $gear
## [1] 4 4 4 3 3 3
## 
## $carb
## [1] 4 4 1 1 2 1

apply() takes as an input dataframe or matrix and applies a function on either rows or columns, returning a single value for each of them. Row-apply and column-apply are specified by MARGIN parameter: 1 for rows and 2 for columns.

apply(mtcars, MARGIN = 1, mean)

##           Mazda RX4       Mazda RX4 Wag          Datsun 710      Hornet 4 Drive   Hornet Sportabout             Valiant          Duster 360           Merc 240D            Merc 230 
##            29.90727            29.98136            23.59818            38.73955            53.66455            35.04909            59.72000            24.63455            27.23364 
##            Merc 280           Merc 280C          Merc 450SE          Merc 450SL         Merc 450SLC  Cadillac Fleetwood Lincoln Continental   Chrysler Imperial            Fiat 128 
##            31.86000            31.78727            46.43091            46.50000            46.35000            66.23273            66.05855            65.97227            19.44091 
##         Honda Civic      Toyota Corolla       Toyota Corona    Dodge Challenger         AMC Javelin          Camaro Z28    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##            17.74227            18.81409            24.88864            47.24091            46.00773            58.75273            57.37955            18.92864            24.77909 
##        Lotus Europa      Ford Pantera L        Ferrari Dino       Maserati Bora          Volvo 142E 
##            24.88027            60.97182            34.50818            63.15545            26.26273

apply(mtcars, MARGIN = 2, mean)

##        mpg        cyl       disp         hp       drat         wt       qsec         vs         am       gear       carb 
##  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750   0.437500   0.406250   3.687500   2.812500

parallel::mclapply(my_vector, mc.cores = 2, FUN = function(x) {
  x ** 2
})

Tidy (long) data

Tidyverse packages are not meant to replace the base R. Especially, they are not suited for matrix manipulation, and are specifically used for long data (also called tidy data). More about tidy data can be found here.

The following three rules make a dataset tidy:

• Each variable must have its own column.
• Each observation must have its own row.
• Each value must have its own cell.

Following three rules makes a dataset tidy: variables are in columns, observations are in rows, and values are in cells.

As you will see, we will be mostly working with biological matrices (= data in wide format). For example, consider this datafame:

(df_wide <- data.frame(sample_1 = 1:2, sample_2 = 3:4, sample_3 = 5:6, row.names = c("gene_A", "gene_B")))

For each sample (column) and each gene (row) we have its value (cell). This form of data is usable for functions operating on matrices (e.g. PCA), but as you will see later, not for e.g. plotting (besides base R plotting).

This dataframe is wide, and thus it is not tidy. You can see that rownames and colnames are, in fact, values of variables Gene and Sample, respectively. So how would this dataframe look like in tidy format? Don’t be surprised now, we will get to this code later.

library(magrittr)
library(tidyverse)

df_long <- df_wide %>%
  tibble::rownames_to_column("Gene") %>%
  tidyr::pivot_longer(-Gene, names_to = "Sample", values_to = "Count")
df_long

Now df_long is tidy and corresponds to long format of the original dataframe. We actually did this:

Converting wide data to long format.

Why is long data useful for some purposes (mainly for plotting)? Well, for biological matrix like the one in df, you usually also have a sample sheet, for example:

(sample_sheet <- data.frame(Sample = c("sample_1", "sample_2", "sample_3"), Age = 1:3, Treatment = c("yes", "no", "yes")))

Would it be neat to have this information alongside the long data (df_long)? I think so:

dplyr::left_join(df_long, sample_sheet, by = "Sample")

As you will see later, this long form of data is extremely useful for plotting. In this blog post, there is a nice example of using tidy data and tidyverse in genomics.

Non-standard evaluation (NSE)

Let’s go ahead a bit and look at this select() function from dplyr package (will be introduced later):

dplyr::select(mtcars, cyl, mpg) %>% head()

You can see that select() is filtering columns of a dataframe. But it has a strange syntax, hasn’t it? We have used variables, which are not defined in our current environment. Don’t be scared by this dark magic - it is based on R’s lazy evaluation (details) and it is called non-standard evaluation (NSE). In short, R only evaluates the expression if the expression is actually used. Just a quick example:

my_fn <- function(x, y, z) {
  print("Hello!")
}

my_fn(a, b, c)

## [1] "Hello!"

You can see we are passing non-existent variables a, b and c to my_fn(), but no errors are thrown! That is because R is evaluating variables (or more generally expressions) only when needed (i.e. lazily).

Also R has quite advanced metaprogramming abilities, so inside a function, we can do things like: “Get the name (symbol) of variable, which was passed to the function.”. So for example, you can find out that you passed variable named a as the argument x and expression mean(b) as the argument y:

my_fn <- function(x, y, z) {
  print(substitute(x))
  print(substitute(y))
}

my_fn(a, mean(b), c)

## a
## mean(b)

Magic! 🧙‍♂️

This is very roughly how most tidyverse functions work. They are using NSE to look first at variables inside a passed object (columns = variables). Most of these functions are designed to have their first parameter as data object. When variables are not found inside the object, function will look for variables in the global environment (this can cause unwanted problems, when object has variable with the same name as the one in global environment). So in the example above, cyl and mpg variables are evaluated in the context of mtcars object and environment defined in select() function call.

Of note, there are some base R functions which also use NSE. Some examples:

cyl2 <- with(mtcars, {cyl ** 2})
head(cyl2)

## [1] 36 36 16 36 64 36

subset(mtcars, hp < 90)

NSE is a quite advanced topic and its knowledge is mainly needed for implementation of custom functions, which share the tidyverse philosophy. See vignette("programming") and friendlyeval, a friendly interface for NSE programming.

Now we will quickly look at the most popular tidyverse packages.

magrittr - a pipe operator

The pipe operator is probably the basic compound of shell programming. Fortunately for us, it is also implemented in R package magrittr.

The pipe operator will become an official part of the R in the next version 4.1.0! Anyways, besides a bit worse performance, magrittr will still offer more flexibility.

The pipe operator (see ?`%>%`) is forwarding output from the left hand side (LHS) to input on the right hand side (RHS). The default behaviour is to pipe output to the first argument of a RHS function.

So instead of using nested functions like

head(rownames(mtcars))

## [1] "Mazda RX4"         "Mazda RX4 Wag"     "Datsun 710"        "Hornet 4 Drive"    "Hornet Sportabout" "Valiant"

you can use the pipe operator:

rownames(mtcars) %>% head()

## [1] "Mazda RX4"         "Mazda RX4 Wag"     "Datsun 710"        "Hornet 4 Drive"    "Hornet Sportabout" "Valiant"

And there are even more types of pipes.

Piped object is available on the RHS as the . object. It is used when you want:

• To use an argument different than first of the RHS function:

c("first", "second", "third") %>% paste("I am", ., sep = " ")

## [1] "I am first"  "I am second" "I am third"

• To use the LHS output as a regular R object:

mtcars %>% .$cyl

##  [1] 6 6 4 6 8 6 8 4 4 6 6 8 8 8 8 8 8 4 4 4 4 8 8 8 8 4 4 4 8 6 8 4

mtcars %>% .[1:5, 1:3]

• To use a lambda expression, which behaves like a function with a single parameter .. Don’t wrap your head around - just wrap an expression on the RHS to {} if you need to use nested functions, e.g. mtcars %>% {head(rownames(.))}

mtcars %>% {head(.$cyl)}

## [1] 6 6 4 6 8 6

mtcars %>% {
  my_data <- .
  c(min(my_data), max(my_data))
}

## [1]   0 472

Or you can define an anonymous function:

mtcars %>% (function(my_data) c(min(my_data), max(my_data)))

## [1]   0 472

Using arithmetic and other operators together with pipe can be more pleasant with provided aliases.

Add and subtract operators are working fine:

1 %>% + 2 %>% - 3

## [1] 0

But for the sake of consistency, rather use operator aliases:

1 %>% add(2) %>% subtract(3)

## [1] 0

It is not possible to use some operators without aliases, e.g. == or *:

1 %>% == 2
1 %>% * 3

## Error: <text>:1:7: unexpected '=='
## 1: 1 %>% ==
##           ^

1 %>% equals(2)

## [1] FALSE

1 %>% multiply_by(3)

## [1] 3

In fact, operators are called infix functions and all aliases are actually shorthands for writing operators as “classic” functions, e.g. `*`(x, y). So the examples above can be also written as:

1 %>% `==`(2)

## [1] FALSE

1 %>% `*`(3)

## [1] 3

But rather don’t do that 😄

tibble - an enhanced data.frame

There are some differences between tibble and base R data.frame (see link above and examples here), but the most visible difference is in printing of a tibble:

tibble::as_tibble(mtcars)

tibble is showing the data types of variables. Also you can see that there aren’t rownames - this is strange, but it is one of the features of tibble. Authors of tibble package are convinced that rownames should be represented by variable (column) and not by object’s attribute, accessible by rownames() method.

This is bad particularly for usage of tidyverse together with Bioconductor, where rownames are an essential component. Many tidyverse functions return a tibble or drop rownames, but most of times we rather want data.frame with defined rownames. tibble can be converted to a data.frame by as.data.frame() function.

To preserve rownames in a new column, use the tibble::rownames_to_column("name_of_new_column") function.

Useful is the tribble() function for more readable creation of (small) tables in code:

my_tibble <- tibble::tribble(
  ~column_A, ~column_B,
  "a",   1,
  "b",   2,
  "c",   3
)
my_tibble

tibble is also useful for handling of nested data.frame-like objects or lists:

my_tibble$column_C <- list(LETTERS, LETTERS, LETTERS)
my_tibble$column_D <- list(mtcars, mtcars, mtcars)
my_tibble

my_tibble$column_C[[1]]

##  [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W" "X" "Y" "Z"

head(my_tibble$column_D[[1]])

dplyr - data manipulation

dplyr is probably the most popular package from tidyverse. It helps to solve the most common data manipulation challenges and its sibling dbplyr also supports remote data stored in a database (MySQL, PostreSQL, SQLite, and others).

Some functions in dplyr have names similar to functions from other packages. Thus, it is recommended to always specify the dplyr:: namespace.

• Cheatsheet

We will use a shortened mtcars dataset for examples:

mtcars_sm <- head(mtcars)

dplyr::select(mtcars_sm, cyl, mpg)

dplyr::select(mtcars_sm, -cyl, -mpg)

dplyr::select(mtcars_sm, disp:drat)

dplyr::select(mtcars_sm, starts_with("c"))

dplyr::select(mtcars_sm, contains("ar"))

dplyr::select(mtcars_sm, last_col(1))

dplyr::select(mtcars_sm, one_of(c("mpg", "cyl", "garbage")))

## Warning: Unknown columns: `garbage`

dplyr::select(mtcars_sm, mpg, cyl, garbage)

## Error: Can't subset columns that don't exist.
## x Column `garbage` doesn't exist.

dplyr::select(mtcars_sm, cylinders = cyl, horse_power = hp)

dplyr::rename(mtcars_sm, cylinders = cyl, horse_power = hp)

dplyr::select(mtcars_sm, cylinders = cyl, horse_power = hp, everything())

dplyr::pull(mtcars_sm, cyl)

## [1] 6 6 4 6 8 6

dplyr::arrange(mtcars_sm, cyl, disp)

dplyr::arrange(mtcars_sm, desc(cyl))
# or
# dplyr::arrange(mtcars_sm, -cyl)

dplyr::filter(mtcars_sm, hp > 100, cyl == 6)

dplyr::filter(mtcars_sm, hp > 100 & cyl == 6)

dplyr::filter(mtcars_sm, hp > 100 | cyl == 6)

dplyr::filter(mtcars_sm, hp > mean(hp) - 10)

tibble::rownames_to_column(mtcars_sm, "car_name") %>%
  dplyr::mutate(cyl2 = cyl * 2) %>%
  dplyr::select(car_name, cyl, cyl2)

dplyr::mutate(mtcars_sm, mpg = NULL, cyl = cyl * 2)

dplyr::mutate(mtcars_sm, car_category = if_else(hp >= 100, "strong", "weak")) %>%
  dplyr::select(hp, car_category)

dplyr::mutate(mtcars_sm, car_category = if_else(hp >= 100, "strong", 50)) %>%
  dplyr::select(hp, car_category)

## Error: Problem with `mutate()` input `car_category`.
## x `false` must be a character vector, not a double vector.
## ℹ Input `car_category` is `if_else(hp >= 100, "strong", 50)`.

dplyr::mutate(
  mtcars,
  car_category = case_when(
    hp < 80 ~ "weak",
    hp >= 80 & hp <= 100 ~ "middle",
    hp > 100 ~ "strong"
  )
) %>%
  dplyr::select(hp, car_category)

dplyr::mutate(mtcars_sm, am_char = recode(am, `0` = "zero", `1` = "one")) %>%
  dplyr::select(am, am_char)

factor(mtcars_sm$am) %>%
  levels()

## [1] "0" "1"

dplyr::mutate(
  mtcars_sm,
  am = factor(am),
  am_char = recode_factor(am, `1` = "one", `0` = "zero")
) %>%
  dplyr::pull(am_char) %>%
  levels()

## [1] "one"  "zero"

v <- c(1, 1, 2, 3, 10)
if_else(v < 2, "< 2", ">= 2")

## [1] "< 2"  "< 2"  ">= 2" ">= 2" ">= 2"

(am_factor <- factor(mtcars_sm$am))

## [1] 1 1 1 0 0 0
## Levels: 0 1

recode(am_factor, `0` = "zero", `1` = "one")

## [1] one  one  one  zero zero zero
## Levels: zero one

recode_factor(am_factor, `1` = "one", `0` = "zero")

## [1] one  one  one  zero zero zero
## Levels: one zero

by_cyl <- dplyr::group_by(mtcars, cyl)
by_cyl

dplyr::summarise(
  by_cyl,
  n = n(),
  mpg_mean = mean(mpg),
  hp_max = max(hp),
  qsec_median = median(qsec),
  wt_sd = sd(wt)
)

dplyr::filter(mtcars, mpg > mean(mpg))

dplyr::filter(by_cyl, mpg > mean(mpg))

mtcars_prices <- tibble::tribble(
  ~cyl, ~price,
  4,    1000,
  6,    5000,
  8,    10000
)

dplyr::left_join(mtcars_sm, mtcars_prices, by = "cyl") %>%
  dplyr::select(cyl, price)

dplyr::bind_rows(mtcars_sm, mtcars_prices) %>%
  dplyr::select(cyl, price, everything())

tidyr - tools for tidy data

The goal of tidyr is to help you to create tidy data. Just to remind you, tidy data is data where:

• Every column is variable.
• Every row is an observation.
• Every cell is a single value.

Cheatsheet on reading and tidying data.

Let’s look at some of basic tasks tidyr can do:

df_wide

df_long <- df_wide %>%
  tibble::rownames_to_column("Gene") %>%
  tidyr::pivot_longer(-Gene, names_to = "Sample", values_to = "Count")
df_long

df_wide2 <- df_wide %>%
  tibble::rownames_to_column("Gene_symbol") %>%
  dplyr::mutate(
    Gene_ID = c("ENSG0001", "ENSG0002"),
    Gene_class = c("protein_coding", "miRNA")
  )

df_wide2

df_long2 <- df_wide2 %>%
  tidyr::pivot_longer(-c(Gene_symbol, Gene_ID, Gene_class), names_to = "Sample", values_to = "Count")

df_long2

df_wide2 %>%
  tidyr::pivot_longer(-starts_with("Gene_"), names_to = "Sample", values_to = "Count")

tidyr::pivot_wider(df_long, names_from = "Sample", values_from = "Count")

tidyr::pivot_wider(df_long2, names_from = "Sample", values_from = "Count")

(df <- data.frame(x = c(NA, "a.b", "a.d", "b.c")))

tidyr::separate(df, x, c("A", "B"), sep = "\\.")

(df <- data.frame(x = c(NA, "a-1-x", "a-2-y", "a-3")))

tidyr::extract(df, x, c("A", "B", "C"), regex = "([a-z])-(\\d)-?([a-z])?")

(df <- data.frame(x = c(NA, "a", "b", "c"), y = c("a", "f", "g", "h")))

tidyr::unite(df, "x_y", x, y, sep = "_", remove = FALSE)

(df <- data.frame(x = c(1, NA, 2), y = c(1, 2, NA)))

tidyr::drop_na(df)

tidyr::drop_na(df, x)

tidyr::replace_na(df, list(x = 5, y = 10))

tidyr::replace_na(c(1, NA, 3), 2)

## [1] 1 2 3

stringr - consistent wrappers for common string operations

stringr is a nice replacement for base R collection of string manipulation functions (e.g. grep(), gsub(), etc., see ?grep). In my opinion, stringr API is more user-friendly. All its functions start with str_ prefix and are vectorised over input string and pattern.

Useful links

• Cheatsheet
• RegExplain - RStudio addin providing a friendly interface for working with regular expressions and functions from stringr.
• Regex101 - cool web application for regex testing.

Because stringr homepage offers a nice collection of examples, I will just copy-paste them here and possibly add more information.

library(stringr)

x <- c("why", "video", "cross", "extra", "deal", "authority")
str_length(x)

## [1] 3 5 5 5 4 9

str_c(…, sep = "", collapse = NULL) concatenates a vector of strings.

str_c(x, collapse = ", ")

## [1] "why, video, cross, extra, deal, authority"

str_sub(string, start = 1L, end = -1L) extracts a substring by a position. Positions can be also negative.

str_sub(x, 1, 2)

## [1] "wh" "vi" "cr" "ex" "de" "au"

str_sub(x, -2)

## [1] "hy" "eo" "ss" "ra" "al" "ty"

str_detect(string, pattern, negate = FALSE) tells you if there’s any match to the pattern. Handy in dplyr::filter().

str_detect(x, "[aeiou]")

## [1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE

str_detect(x, "[aeiou]", negate = TRUE)

## [1]  TRUE FALSE FALSE FALSE FALSE FALSE

str_count(x, pattern) counts the number of patterns.

str_count(x, "[aeiou]")

## [1] 0 3 1 2 2 4

str_subset(string, pattern, negate = FALSE) extracts the matching components. str_which(string, pattern, negate = FALSE) gives indices of the matching components.

str_subset(x, "[aeiou]")

## [1] "video"     "cross"     "extra"     "deal"      "authority"

str_subset(x, "[aeiou]", negate = TRUE)

## [1] "why"

str_which(x, "[aeiou]")

## [1] 2 3 4 5 6

str_locate(string, pattern) gives the position of the match.

str_locate(x, "[aeiou]")

##      start end
## [1,]    NA  NA
## [2,]     2   2
## [3,]     3   3
## [4,]     1   1
## [5,]     2   2
## [6,]     1   1

str_extract(string, pattern) extracts the text of the match and returns a vector.

x

## [1] "why"       "video"     "cross"     "extra"     "deal"      "authority"

str_extract(x, "[aeiou]")

## [1] NA  "i" "o" "e" "e" "a"

str_match(string, pattern) extracts parts of the match defined by parentheses (capturing groups). Return value = matrix, with rows being input vector values, and columns being capturing groups. First column represents the full match.

x

## [1] "why"       "video"     "cross"     "extra"     "deal"      "authority"

str_match(x, "(.)[aeiou](.)")

##      [,1]  [,2] [,3]
## [1,] NA    NA   NA  
## [2,] "vid" "v"  "d" 
## [3,] "ros" "r"  "s" 
## [4,] NA    NA   NA  
## [5,] "dea" "d"  "a" 
## [6,] "aut" "a"  "t"

str_replace(string, pattern, replacement) replaces the matches with new text. replacement also supports references of the form \1, \2, etc., which will be replaced with the contents of the respective matched group (created by ()). str_remove(string, pattern) is an alias for str_replace(string, pattern, "").

str_replace(x, "[aeiou]", "?")

## [1] "why"       "v?deo"     "cr?ss"     "?xtra"     "d?al"      "?uthority"

str_replace_all(x, "[aeiou]", "?")

## [1] "why"       "v?d??"     "cr?ss"     "?xtr?"     "d??l"      "??th?r?ty"

str_replace(x, "(a)", "\\1_")

## [1] "why"        "video"      "cross"      "extra_"     "dea_l"      "a_uthority"

str_remove(x, "[aeiou]")

## [1] "why"      "vdeo"     "crss"     "xtra"     "dal"      "uthority"

str_remove_all(x, "[aeiou]")

## [1] "why"   "vd"    "crss"  "xtr"   "dl"    "thrty"

str_split(string, pattern, n = Inf, simplify = FALSE) splits up a string into multiple pieces.

str_split(c("a,b", "c,d,e"), ",")

## [[1]]
## [1] "a" "b"
## 
## [[2]]
## [1] "c" "d" "e"

glue - string interpolation

glue offers the same functionality, which you already know from e.g. Bash ("$variable") or Python (f"{variable}").

library(glue)
name <- "Fred"
age <- 50
anniversary <- as.Date("1991-10-12")
glue("My name is {name}, my age next year is {age + 1}, my anniversary is {format(anniversary, '%A, %B %d, %Y')}.")

## My name is Fred, my age next year is 51, my anniversary is Saturday, October 12, 1991.

Equivalently, we can break long string to multiple parameters on separate lines and use named arguments to assign temporary variables:

glue(
  "My name is {name},",
  " my age next year is {age + 1},",
  " my anniversary is {format(anniversary, '%A, %B %d, %Y')}.",
  name = "Joe",
  age = 40,
  anniversary = as.Date("2001-10-12")
)

purrr - functional programming tools

purrr is basically a handy replacement for lapply() and sapply(), but if you go in more depth, it offers much more functionality. We will just look at the most basic functions.

• Cheatsheet

map() is similar to lapply(), but allows to use a special formula notation for functions. In this formula, iterated variable is referenced by . (dot). In the example below, ~ . + 1 is equivalent to function(x) {x + 1}.

my_list <- list(first = 1:3, second = 1:10, third = 20)
my_list

## $first
## [1] 1 2 3
## 
## $second
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $third
## [1] 20

purrr::map(my_list, ~ . + 1)

## $first
## [1] 2 3 4
## 
## $second
##  [1]  2  3  4  5  6  7  8  9 10 11
## 
## $third
## [1] 21

map() or pluck() can be used to retrieve a specific subelement of each element. Here, we return the first element of each nested vector:

purrr::map(my_list, 1)

## $first
## [1] 1
## 
## $second
## [1] 1
## 
## $third
## [1] 20

Similarly, if nested lists or vectors are named, we can use their names instead of integer position.

my_named_list <- list(
  list(first = 1:2, second = "Hello", third = 20),
  list(first = 5, second = 1)
)
my_named_list

## [[1]]
## [[1]]$first
## [1] 1 2
## 
## [[1]]$second
## [1] "Hello"
## 
## [[1]]$third
## [1] 20
## 
## 
## [[2]]
## [[2]]$first
## [1] 5
## 
## [[2]]$second
## [1] 1

purrr::map(my_named_list, "second")

## [[1]]
## [1] "Hello"
## 
## [[2]]
## [1] 1

map() always returns a list, while its map_<type>() versions returns a vector of specified type (e.g. character - map_chr()).

purrr::map_dbl(my_list, mean)

##  first second  third 
##    2.0    5.5   20.0

Using a predicate function, keep() and discard() will filter elements:

purrr::keep(my_list, ~ all(. < 5))

## $first
## [1] 1 2 3

purrr::discard(my_list, ~ all(. < 5))

## $second
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $third
## [1] 20

purrr::keep(my_named_list, ~ length(.) > 2)

## [[1]]
## [[1]]$first
## [1] 1 2
## 
## [[1]]$second
## [1] "Hello"
## 
## [[1]]$third
## [1] 20

Using dynamic variable names in aesthetics

You probably wondered, how can you use dynamic variable names in aes(). For example, consider a function plot_pca(data, color_by) which calculates PCA and shows scatterplot colorized by a specified variable. Although it is possible via NSE (see vignette("programming")), the easiest way is to use aes_string() function.

color_by <- "cyl"

ggplot(mtcars) +
  geom_point(aes_string(x = "mpg", y = "qsec", color = color_by))

The disadvantage is that you cannot use the capabilities of NSE. In the example above, you cannot specify factor(color_by) to treat cyl as a factor variable.

Packages extending the ggplot2

library(patchwork)
p1 <- ggplot(mtcars) + geom_point(aes(mpg, disp))
p2 <- ggplot(mtcars) + geom_boxplot(aes(gear, disp, group = gear))
p3 <- ggplot(mtcars) + geom_smooth(aes(disp, qsec))
p4 <- ggplot(mtcars) + geom_bar(aes(carb))

p1 + p2

p1 / p2

p1 + p2 + p3 + p4

p1 + p2 + p3 + p4 + plot_layout(nrow = 3, byrow = FALSE)

p1 + (p2 / p3)

p1 / (p2 + p3)

(p1 | (p2 / p3)) + 
  plot_annotation(title = "The surprising story about mtcars", tag_levels = "A")

# Or you can pass a list of plots:
# plots <- list(p1, p2, p3, p4) 
# wrap_plots(plots)
wrap_plots(p1, p2, p3, p4)

wrap_plots(p1, p2, p3, p4, nrow = 1)

ggpubr::ggboxplot(
  mtcars, x = "cyl", y = "qsec",
  color = "cyl", shape = "cyl",
  add = "jitter"
)

GGally::ggpairs(mtcars, columns = c("cyl", "mpg", "qsec"))

ggplot(mtcars) +
  geom_point(aes(x = hp, y = qsec, color = factor(gear))) +
  ggforce::facet_zoom(x = gear == 4) +
  theme_bw()

corr <- round(cor(mtcars), 1)
corr[1:5, 1:5]

##       mpg  cyl disp   hp drat
## mpg   1.0 -0.9 -0.8 -0.8  0.7
## cyl  -0.9  1.0  0.9  0.8 -0.7
## disp -0.8  0.9  1.0  0.8 -0.7
## hp   -0.8  0.8  0.8  1.0 -0.4
## drat  0.7 -0.7 -0.7 -0.4  1.0

p_mat <- ggcorrplot::cor_pmat(mtcars)
p_mat[1:5, 1:5]

##               mpg          cyl         disp           hp         drat
## mpg  0.000000e+00 6.112687e-10 9.380327e-10 1.787835e-07 1.776240e-05
## cyl  6.112687e-10 0.000000e+00 1.802838e-12 3.477861e-09 8.244636e-06
## disp 9.380327e-10 1.802838e-12 0.000000e+00 7.142679e-08 5.282022e-06
## hp   1.787835e-07 3.477861e-09 7.142679e-08 0.000000e+00 9.988772e-03
## drat 1.776240e-05 8.244636e-06 5.282022e-06 9.988772e-03 0.000000e+00

p_corr <- ggcorrplot::ggcorrplot(corr, method = "circle", p.mat = p_mat)
print(p_corr)

ggplot(mtcars2, aes(x = qsec, y = mpg)) +
  geom_point(aes(color = factor(gear), size = hp)) +
  ggrepel::geom_text_repel(aes(label = rownames(mtcars)), color = "black")

Additional themes

p1 <- ggplot(mtcars) +
  geom_point(aes(x = hp, y = qsec, color = factor(cyl))) +
  theme_bw() +
  ggsci::scale_color_npg()

p2 <- ggplot(mtcars) +
  geom_histogram(aes(x = qsec, fill = factor(cyl)), colour = "black", binwidth = 1, position = "dodge") +
  theme_bw() +
  ggsci::scale_fill_npg()

p1 + p2

ggplot(mtcars) +
  geom_point(aes(x = hp, y = qsec, color = factor(cyl))) +
  hrbrthemes::theme_modern_rc()

ggplot(mtcars) +
  geom_point(aes(x = hp, y = qsec, color = factor(cyl))) +
  see::theme_abyss()

janitor - table summaries, duplicates etc.

Personally, I am mainly using these functions:

tabyl() (vignette) creates a frequency table. It is a tidy replacement for base R table():

table(mtcars$cyl)

## 
##  4  6  8 
## 11  7 14

janitor::tabyl(mtcars, cyl)

table(mtcars$cyl, mtcars$gear)

##    
##      3  4  5
##   4  1  8  2
##   6  2  4  1
##   8 12  0  2

janitor::tabyl(mtcars, cyl, gear)

table(mtcars$cyl, mtcars$gear, mtcars$am)

## , ,  = 0
## 
##    
##      3  4  5
##   4  1  2  0
##   6  2  2  0
##   8 12  0  0
## 
## , ,  = 1
## 
##    
##      3  4  5
##   4  0  6  2
##   6  0  2  1
##   8  0  0  2

janitor::tabyl(mtcars, cyl, gear, am)

## $`0`
##  cyl  3 4 5
##    4  1 2 0
##    6  2 2 0
##    8 12 0 0
## 
## $`1`
##  cyl 3 4 5
##    4 0 6 2
##    6 0 2 1
##    8 0 0 2

adorn_*() append row and/or column with summary statistics:

janitor::tabyl(mtcars, cyl, gear) %>%
  janitor::adorn_totals(c("row", "col")) %>%
  janitor::adorn_percentages("row") %>%
  janitor::adorn_pct_formatting(rounding = "half up", digits = 0) %>%
  janitor::adorn_ns() %>%
  as.data.frame()

get_dupes() returns duplicated rows:

janitor::get_dupes(mtcars, mpg, hp)

plotly - interactive HTML plots

plotly is a JavaScript library with wrappers for R and Python. It is somewhat similar to ggplot2 in terms of The Grammar of Graphics. Output is a HTML file or in case of RMarkdown, it can be directly embedded as a chunk output.

We don’t go into details, so if you are interested, see examples. We will try to replicate the mtcars scatterplot done in ggplot2:

p_ggplot2 <- ggplot(mtcars2, aes(x = qsec, y = mpg)) +
  geom_point(aes(color = factor(gear), size = hp)) +
  theme_bw()

print(p_ggplot2)

p_plotly <- plotly::plot_ly(
  data = mtcars,
  x = ~qsec, y = ~mpg, color = ~factor(gear), size = ~hp,
  type = "scatter", mode = "markers",
  text = ~glue("mpg: {mpg}<br>qsec: {qsec}<br>gear: {gear}<br>hp: {hp}")
)

p_plotly

plotly can also convert ggplot2 objects, but it is not always accurate:

plotly::ggplotly(p_ggplot2)

But this seems to be quite good, even better than p_plotly ☺️

plotly produces objects inheriting from htmlwidget class. We can use htmlwidgets package to save plotly object into HTML file:

class(p_plotly)

## [1] "plotly"     "htmlwidget"

htmlwidgets::saveWidget(p_plotly, "p_plotly.html")

p_plotly.html

heatmaply - interactive HTML heatmaps

heatmaply uses plotly underneath. Its API is mostly compatible with the classic heatmap.2() function.

heatmaply::heatmaply(mtcars)

It has also functionality similar to ggcorrplot:

cor <- psych::corr.test(mtcars)

heatmaply::heatmaply_cor(
  cor$r,
  node_type = "scatter",
  point_size_mat = -log10(cor$p), 
  point_size_name = "-log10(p-value)",
  label_names = c("x", "y", "Correlation")
)

ComplexHeatmap

ComplexHeatmap is probably the most complex package for heatmap visualization, while still providing very user-friendly approach.

• Reference book
• Paper: Complex heatmaps reveal patterns and correlations in multidimensional genomic data

ComplexHeatmap was initially created for complex multi-omics datasets, combining e.g. gene expression with methylation levels. Its design follows the object-oriented scheme, so heatmap body, annotations, legends etc. are self-contained building blocks. It can do a lot and it can be a bit overhelming, but for simple heatmaps (e.g. z-score of logcounts), unlike heatmap.2 and similar packages, ComplexHeatmap is more user-friendly and provides you a great assistance.

Some time ago, I created an introduction for my labmates, which is available in this GitHub repository, but you can find all the necessary files in E02-intro_to_R/ComplexHeatmap-intro: there is a source Rmd file (ComplexHeatmap.Rmd) and its rendered HTML document (ComplexHeatmap.html).

Although there are many popular heatmap packages, I will deliberately not mention them here, as they have many drawbacks. In my opinion, ComplexHeatmap is currently the best heatmap package, especially for genomics and transcriptomics.

BiocParallel

BiocParallel package provides easy-to-use functions for parallel execution. We won’t go into details and just introduce probably the most used function bplapply(), a parallelized version of lapply().

As you know, lapply(X, FUN, ...) applies FUN with optional arguments ... to each element of X:

x <- 1:5
fn <- function(i) i^2
lapply(x, fn)

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 4
## 
## [[3]]
## [1] 9
## 
## [[4]]
## [1] 16
## 
## [[5]]
## [1] 25

bplapply(X, FUN, ..., BPREDO = list(), BPPARAM = bpparam()) does the same, but elements are split to chunks, which are processed in parallel. Result is guaranteed to have the same order of elements as in X.

BiocParallel::bplapply(x, fn)

## [[1]]
## [1] 1
## 
## [[2]]
## [1] 4
## 
## [[3]]
## [1] 9
## 
## [[4]]
## [1] 16
## 
## [[5]]
## [1] 25

By default, all available cores are used. To control for resources, use the BPPARAM argument:

(bpparam <- BiocParallel::MulticoreParam(workers = 2))

## class: MulticoreParam
##   bpisup: FALSE; bpnworkers: 2; bptasks: 0; bpjobname: BPJOB
##   bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
##   bpRNGseed: ; bptimeout: 2592000; bpprogressbar: FALSE
##   bpexportglobals: TRUE
##   bplogdir: NA
##   bpresultdir: NA
##   cluster type: FORK

BiocParallel::bplapply(x, fn, BPPARAM = bpparam)

Many functions in BioConductor packages will offer you an argument to pass a BPPARAM object. If you are working on some high-performance server, this can speed up computations a lot.

Parallel processing is hard to debug, even with try-catch equivalent function bptry(). Also, you can’t use browser() inside FUN. But sometimes is just enough to run bplapply() in non-parallel (serial) mode. This is possible with SerialParam().

(bpparam <- BiocParallel::SerialParam())
fn <- function(i) {
  browser()
  i^2
}

BiocParallel::bplapply(x, fn, BPPARAM = bpparam)

BiocParallel supports jobs in distributed systems and much more, but for now we can live only with the simple multicore usage.